#############################################
#############################################
### REPLICATION MATERIAL                  ###  
### Gender Quotas and Political Knowledge ###
### Quotas Knowledge Analysis (main)      ###
#############################################
#############################################

# Libraries needed
library(tidyverse)
library(stargazer)
library(haven)
library(lme4)
library(rstanarm)
library(tidybayes)
library(BayesPostEst)
library(ggplot2)
library(ggridges)
library(magrittr)
library(dplyr)
library(purrr)
library(forcats)
library(tidyr)
library(modelr)
library(ggdist)
library(tidybayes)
library(cowplot)
library(rstan)
library(RColorBrewer)
library(gganimate)
library(wrangle)

# Set plot theme
theme_set(theme_tidybayes() + panel_border())

# Set seed and stan options (to make it run faster)
set.seed(123)
options(mc.cores = parallel::detectCores())

# Read RDS data
analysis.2 <- readRDS("data/analysis/analysis_2_df.rds")

# Modify full dataset
analysis.int <- analysis.2 %>%
  mutate(female = recode(gender,
                         "1" = 0,
                         "0" = 1)) %>%
  select(ID, country, name, year, know_q_num, know_correct, 
         female, quota, age, education, GDP_growth, counter) %>%
  mutate(age = cut(age, 5,labels = FALSE)) %>%
  mutate(education = cut(education, 5, labels = FALSE)) %>%
  mutate(treated = ifelse(name == "Belgium" | name == "France" | name == "Greece" |
                            name == "Ireland" | name == "Portugal" |
                            name == "Spain" , 1, 0)) %>%
  write_rds("data/analysis/analysis_int.rds") %>%
  glimpse()


# Binomial version of the dataset - full sample
analysis.int.binomial <- analysis.int %>%
  group_by(female, quota, age, treated, counter, education, 
           GDP_growth, country, year, know_q_num) %>%
  summarize(n_correct = sum(know_correct),
            n_incorrect = n() - sum(know_correct)) %>%
  write_rds("data/analysis/analysis_int_binomial.rds") %>%
  glimpse()

###################################
###################################
## MODEL AND RESULTS TABLE (A2) ###
###################################
###################################

fit.1 <- stan_glmer(cbind(n_correct, n_incorrect) ~ female*quota*age +
                        education + 
                        GDP_growth + 
                        (1 | country) + 
                        (1 | year) + 
                        (1 | know_q_num) + 
                        (1 + female | country:know_q_num), 
                      data = analysis.int.binomial,
                      family=binomial(link="logit"),
                      iter = 6000,
                      control = list(adapt_delta = 0.99))

#######################
### Summary - TABLE A2
#######################

print(fit.1, digits = 2)
mcmcReg(fit.1, format = "latex")
mcmcTab(fit.1)

#######################
### Convergence checks
#######################

# Rhat - convergence ok
dev.copy(png,'rhat_fit_1.png')
plot(fit.1, "rhat")
dev.off()

# Prior summary
prior_summary(fit.1)


#################################
### PREDICTED PROBS DATASETS  ###
#################################

## Dataset with predicted probabilities (no random effects)
pred_fit_1_noRE <- crossing(country = unique(analysis.int.binomial$country), 
                              quota = c(0,1), 
                              female = c(0,1),
                              know_q_num = c(1:10),
                              age = c(1:5),
                              education = mean(analysis.int.binomial$education, 
                                               na.rm = T),
                              GDP_growth = mean(analysis.int.binomial$GDP_growth, 
                                                na.rm = T)) %>%
  add_linpred_draws(fit.1, 
                    scale = "response", 
                    re_formula = NA) %>% # analogous to predict()
  write_rds("data/analysis/pred_fit_1_noRE.rds")

## Dataframe with draws
draws <- fit.1 %>%
  spread_draws(b[term,group]) %>%
  filter(str_detect(group, "country:know_q_num:")) %>%
  separate(group, c("country_2", "q_num_2", "country", "q_num"), ":") %>%
  write_rds("/data/analysis/draws_fit_1.rds")

#Dataframe with predicted probabilities (with random effects)
pred_fit_1_c_RE <- crossing(country = unique(analysis.int.binomial$country), 
                              quota = c(0,1), 
                              female = c(0,1),
                              year = 2005,
                              know_q_num = c(1:10),
                              age = c(1:5),
                              education = mean(analysis.int.binomial$education, 
                                               na.rm = T),
                              GDP_growth = mean(analysis.int.binomial$GDP_growth, 
                                                na.rm = T)) %>%
  add_linpred_draws(fit.1, 
                    scale = "response", 
                    re_formula = ~ (1 + female | country:know_q_num)) %>% 
  write_rds("/data/analysis/pred_fit_1_c_RE.rds")

#############
#############
### PLOTS ###
#############
#############

##################################################
### Figure 1: Predicted Probabilities - Quotas ###
################################################## 

plot_1 <- readRDS("data/analysis/pred_fit4_noRE.rds")

plot_1$female <- factor(plot_1$female,
                        labels = c("Male", "Female"))
plot_1$quota <- factor(plot_1$quota, 
                       labels = c("No Quotas", "Quotas"))
plot_1$age <- factor(plot_1$age, 
                     labels = c("15-31", "32-48", "49-65", "66-82", "83-99"))

plot_1 %>%
  ggplot(aes(x = female, y = .value, fill = female, color = female)) +
  stat_halfeye() +
  facet_grid(quota ~ .) +
  coord_flip() +
  scale_color_grey(name = "Gender",
                   labels = c("Male", "Female")) +
  scale_fill_grey(name = "Gender",
                  labels = c("Male", "Female")) +
  labs(x = "Gender", y = "Predicted Probabilities")
ggsave("plots/main/fig-1-greyscale.png")


#########################################################
### Figure 2: Predicted Probabilites - Age and quotas ###
######################################################### 

plot_2 <- readRDS("data/analysis/pred_fit4_noRE.rds")

plot_2$female <- factor(plot_2$female,
                        labels = c("Male", "Female"))
plot_2$quota <- factor(plot_2$quota, 
                       labels = c("No Quotas", "Quotas"))
plot_2$age <- factor(plot_2$age, 
                     labels = c("15-31", "32-48", "49-65", "66-82", "83-99"))

plot_2 %>%
  ggplot(aes(x = age, y = .value, fill = female, color = female)) +
  stat_halfeye() +
  facet_grid(quota ~ .) +
  coord_flip() +
  scale_color_grey(name = "Gender",
                   labels = c("Male", "Female")) +
  scale_fill_grey(name = "Gender",
                  labels = c("Male", "Female")) +
  labs(x = "Age Groups", y = "Predicted Probabilities")
ggsave("plots/main/fig-2-greyscale.png")

#########################################################
### Figure 3: Percent Change- Gender Gap - by age     ###
######################################################### 

plot_3 <- readRDS("data/analysis/pred_fit5_c_RE.rds")

percent_change <- plot_3 %>%
  ungroup() %>% 
  mutate(quota = case_when(quota == 1 ~ "quota",
                           quota == 0 ~ "noquota"),
         female = case_when(female == 0 ~ "man",
                            female == 1 ~ "woman")) %>%
  select(-.row) %>% 
  pivot_wider(names_from = c(quota, female), values_from = .value) %>%
  mutate(second_difference = (quota_woman - quota_man) - (noquota_woman - noquota_man)) %>%
  mutate(difference_noquota = (noquota_woman - noquota_man)) %>%
  mutate(percent_change = second_difference/difference_noquota) %>%
  mutate(quotas = ifelse(country == 1 | country == 2 | country == 8 | country == 10 | country == 11 | country == 12, 1, 0)) %>%
  glimpse()

percent_change$country <- factor(percent_change$country,
                                 labels = c("France", "Belgium", 
                                            "Netherlands", "Germany", "Italy",
                                            "Luxembourg", "Denmark", "Ireland",
                                            "United Kingdom", "Greece", "Spain",
                                            "Portugal"))

percent_change$age <- factor(percent_change$age, 
                             labels = c("15-31", "32-48", "49-65", "66-82", "83-99"))

percent_change %>% 
  ggplot(aes(x = percent_change, 
             y = reorder(age, percent_change),
             color = age)) + 
  stat_pointinterval() +
  coord_cartesian(xlim=c(-0.6,0.6)) +
  scale_color_grey(name = "Age Groups",
                   labels = c("15-31", "32-48", "49-65", "66-82", "83-99")) +
  geom_vline(xintercept = 0) +
  labs(x = "Percent Change",
       y = "Age")
ggsave("plots/main/fig-3-grayscale.png")



